function bspline1

%  Solves the BVP below using B-splines 
%       y'' + p(x)y' + q(x)y= f(x)   for xL < x < xr  '
%  where 
%      y(xl) = yL  and y(xR) = yR

%  p=1, q=-2, f=-3*exp(-2*x)
%  xL=0, yL=0  and  xR=1, yR=exp(-2)

% clear all previous variables and plots
clear *
clf

% set boundary conditions
	xL=0; yL=0;
	xR=1; yR=exp(-2);

% calculate and plot exact solution
xe=linspace(xL,xR,100);
for k=1:100
	ye(k)=xe(k)*exp(-2*xe(k));
end;
plot(xe,ye,'-k')
hold on

% define title and axes used in plot
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
title('B-Splines vs Exact Solution','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
box on
axis([0 1 0.0 0.22])
set(gca,'ytick',[0 0.1 0.2]);
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 

% plot numerical solution for various N values
for in=1:2
	N=3*in;

	% generate the points along the x-axis, x(1)=xL and x(n)=xR
	x=linspace(xL,xR,N);
	h=x(2)-x(1);
	
	% calculate the coefficients for B-splines
	a=zeros(1,N); b=zeros(1,N); c=zeros(1,N); z=zeros(1,N);
	hh=h*h;
	for i=1:N
		a(i)=4*(-3+hh*q(x(i)));
		b(i)=6-3*h*p(x(i))+hh*q(x(i));
		c(i)=6+3*h*p(x(i))+hh*q(x(i));
		z(i)=6*hh*f(x(i));
	end;
	z(1)=z(1)-6*yL*b(1); a(1)=a(1)-4*b(1);  c(1)=c(1)-b(1);
	z(N)=z(N)-6*yR*c(N); a(N)=a(N)-4*c(N);  b(N)=b(N)-c(N);
	
	% solve the tri-diagonal matrix problem
	w=tri(a,b,c,z); 
	wL=6*yL-4*w(1)-w(2); wR=6*yR-w(N-1)-4*w(N);
	ww=[wL, w, wR];

	% plot solution using B-spline functions
	np=100;
	xp=linspace(xL,xR,np);
	for ix=1:np
		sum=0;
		for k=1:N+2
			xk=-h+h*(k-1);
			xx=(xp(ix)-xk)/h;
			sum=sum+ww(k)*bspline(xx);
		end;
		yp(ix)=sum;
	end;

	if in==1
		plot(xp,yp,'--r','LineWidth',1)
		legend(' Exact',' N = 3',4)
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
		pause
	elseif in==2
		plot(xp,yp,'-.b','LineWidth',1)
		legend(' Exact',' N = 3',' N = 6',4)
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 
	end;
end
hold off


function g=q(x)
g=-2;

function g=p(x)
g=1;

function g=f(x)
g=-3*exp(-2*x);


function y=bspline(x)
% Calculate the value of a cubic B-spline at point x
x=abs(x) ;
if x > 2,
  y=0 ;
else
  if x > 1,
    y=(2-x)^3/6 ;
  else
    y=2/3-x^2*(1-x/2) ;
  end ;
end ;

% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end